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Abstract — The aim of this paper is to review some of the models 
and solution techniques used in the simulation of high-power 
semiconductor lasers and to address open questions. We discuss 
some of the peculiarities in the description of the optical field 
of wide-aperture lasers. As an example, the role of the substrate 
as a competing waveguide in GaAs-based lasers is studied. The 
governing equations for the investigation of modal instabilities 
and filamentation effects are presented and the impact of the 
thermal-lensing effect on the spatiotemporal behavior of the 
optical field is demonstrated. We reveal the factors that limit the 
output power at very high injecton currents based on a numerical 
solution of the thermodynamic based drift-diffusion equations 
and elucidate the role of longitudinal spatial holeburning. 



I. Introduction 

TREMENDOUS advancements have been achieved in the 
development of high-power semiconductor lasers during 
the last two decades, due to improved cavity design, crystal 
growth, facet passivation, and device cooling technologies 
|[T]-||6). The output power has increased by one order of 
magnitude, the spectral linewidth has decreased thanks to 
Bragg gratings integrated into the cavity, and the beam quality 
has greatly improved thanks to tapered gain regions. Semi- 
conductor lasers are under the way not longer to serve only 
as optical pumps for solid-state lasers but to replace them 
because of their ease-of-use, compactness, efficiency and high 
reliability. However, many of the parameters of high-power 
semiconductor lasers are still far from what is achievable from 
solid-state lasers. 

For example, the reliable maximum output power is cur- 
rently limited to values below 20 W from a 100-^m stripe 
width device and a further increase remains an open challenge. 
High-power semiconductor lasers are also plagued by the 
appearance of multiple peaks in the lateral far- and near field 
profiles as well as in the optical spectrum in particular at 
high injection currents. These effects are caused by modal 
instabilities or by lasing filaments formed by a self-focusing 
mechanism where the refractive index locally increases in 
regions of high optical intensity. 

In order to assess the root causes of these limitations 
physics-based modeling and numerical simulation is more 
important than ever The aim of this paper is to review 
some of the models and solution techniques used in the 
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simulation of high-power semiconductor lasers and to address 
open questions. 

The paper is organized as follows. We will start in Section 
II with a model for the optical field and introduce concepts 
such as the paraxial parabolic equation, roundtrip operator and 
modes. In Section III we discuss the modeling of the nonlinear 
interaction between the optical field and the injected carriers in 
high-power lasers leading to a non-stationary behavior of the 
field pattern. The simulation of the stationary electro-optical 
characteristics of high-power lasers using the thermodynamic 
based drift-diffusion (or energy transport) model will be pre- 
sented in Section IV. We end with a summary and outlook. 

II. Description of the optical field 

A. Paraxial parabolic equation 

The optical field in edge-emitting semiconductor lasers can 
be well described by the paraxial parabolic equations ||7] 

ikonn.dE^ .dE^ 1 
2/3o 



c/3o dt ~^ dz 



E"^ ^0 (1) 



2/3o 2/3o 

for the right and left traveling waves £'+ and E^, respectively. 
Eq. ([T]i can be derived from Maxwell equations with the 



Ansatz 



Eirt,z,t) = e,[E+{ruz,t)e-'^°' 



+ E {rt,z,t)e 



(2) 



for the transverse electric field in order to remove the rapid 
variations with respect to the time t and the coordinate z along 
the cavity axis, and a corresponding Ansatz for the complex 
dielectric function. 



e{rt,z,t,uj)\^g = n^{rt,z,t) 



(3) 



In (T]|3i, Tf — {x,y) are the transverse coordinates. At is the 
transverse Laplacian, cjq and /?o are the real-valued reference 
frequency and reference propagation factor, respectively, fco = 
ujq/c = 27r/Ao with Ao being the reference vacuum wave- 
length and c the vacuum velocity of light, rig = d{ujn)/duj\^g 
is the group index and are the Fourier coefficients of 
a Bragg grating integrated into the cavity. We have further 
assumed a TE-polarized field with e.j., being the unity vector 
and that the dielectric function varies slowly with respect to 
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X. The distribution of the complex-valued refractive index n 
(which may include a time dependence much slower than the 
variation given by I/wq) can be written as 



■ 9 -a 
' 2fco 



(4) 



with g being the gain due to inter band transition and a the 
absorption coefficient due to intraband transitions (e.g. free 
carrier absorption). 

The temporal dispersion of the dielectric function has been 
taken into account only up to the first order in the real part 
of the refractive index n, higher order dispersion and the 
dispersion of the imaginary part have to be properly added 
(see Section pHi. We should note that in a numerical evaluation 
of ([TJ the factor in front of the time-derivative must be 
approximated by 

c/3o Vg 

with a real- valued group velocity Vg. 

Eqs. ([TJ must be supplemented by appropriate boundary 
conditions. At the plane facets of the laser 



E+{0)-roE-iO)=0, 







(6) 



hold at z = and z = L, respectively. We should mention that 
in the paraxial approximation the facet reflectivities rg and 
are input parameters which have to be calculated in advance. 

At the transverse boundary denoted by F one can assume, 
for example, decaying fields or a perfect electric wall. 



lim E"^ =0 or E'^lr.er 



0, 



(7) 



respectively. If only a part of the cross-section of the cavity is 
simulated, a non-reflecting or transparent boundary condition 
has to be used which models the fact that only outgoing waves 
should be present. This boundary condition can be formulated 
only in operator form for the general case. 



dn 



\rter 



-iBE^ 



\rter 



(8) 



where D is the operator for a so-called Dirichlet-to-Neumann 
map |8 1. We will consider below in Subsection II-F an example 
for a one-dimensional case. For spatially and temporally 
constant n, D can be obtained by a factorization of ([TJ |9|. 
A very popular method to implement a boundary condition of 
type ([8]) is the introduction of a so-called perfectly matched 
layer (PML) |10|. However, in the frequency domain it results 
in a large number of spurious modes fTT) . 



subject to the boundary conditions (|6]l-([8]l. The non-trivial so- 
lutions of (j9| may dependent on time via the time-dependence 
of n. The complex-valued relative mode frequencies il„i 
the eigenvalues and the mode profiles E^{rt, z,t) are the 
eigenfunctions of (j9]l. Some mathematical properties and their 
physical impact will be discussed below. The real parts of ftm 
give the wavelengths relative to the reference wavelength Aq, 



dX 

did 



(10) 



and the imaginary parts describe the damping of the modes. 
For a passive cavity, lm{fl„i) > must hold. Lasing modes 
of an active cavity are distinguished by vanishing damping, 
Im(ri„j) = 0, due to the balance of the outcoupling and 
internal losses and the gain. 

The cavity modes fulfill the orthogonality relation 



J nug [E+E„^, + E„^E+,] dxdydz = 



for m ^ m! 

(11) 

which is proven in the Appendix. It is similar to corresponding 
orthogonality relations for cavity modes derived in [12] and 
| fT3| from Maxwell and Helmholtz equations, respectively. 
Note that the integral 



11 



does not define a scalar product 
p4|. In fact, (|9|) defines a non-Hermitian eigenvalue problem, 
because the frequencies Q,m are complex-valued. It is well- 
known from the theory of non-Hermitian operators, that in de- 
pendence of some parameter(s) so-called "exceptional points" 
exist where the eigenvalues (both real and imaginary parts) 
cross and the eigenfunctions become identical |15|. Mode de- 
generacies related to exceptional points have been discovered 
for unstable laser cavities, cf. p6[ and the references therein, 
multi-section lasers f\n\, fTSl and complex planar waveguides 
y_9J. Due to the mode degeneracy, at an exceptional point the 
system of eigenfunctions is no longer complete, but a system 
including the generalized eigenfunctions is, cf. |20|. 

At the exceptional point, the modes remain orthogonal in 
the sense of (111, so that the integral J nngE:^E:^dxdydz 
vanishes. As a consequence, Petermann's K factor pT[ 



K„ 



Jnug [\E+\^ + \E-\^] dxdydz 



2 J nUgEmEmdxdydz 



(12) 



approaches infinity. It is known, that the K factor causes an 
enhancement of the spontaneous emission, which results in a 
broadening of the spectrum of a multi-longitudinal-mode laser 
and the spectral linewidth of a single mode as summarized in 



|22|. Its time-dependence is needed to explain the dynamical 
behavior of multi-section lasers 1231. 



B. Cavity modes 

The cavity modes are time-periodic solutions of the form 
E^e.x^{iQ,jnt) and obey the equations 



konUgVl^ . ,dE^ 



1 



c/3o 



dz 2/3o 

2^2 «2 



2/3o 



^tE^ 



2f± 



2/3o 



-E^^Q (9) 



C. Beam propagation method and roundtrip operator 

For a solution of (|9]l it is convenient to introduce the 
operator 



H(n,z) 
with the abbreviation 

A/3 = 



1 

Wo 



At + A/3 



PI 



2/3o 



(13) 



(14) 
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If f * =0 holds (Fabry-Perot cavity), and if the index n Nowadays, there are more stable and efficient methods to solve 



depends only on the transverse coordinates, 
solution of (j9]l can be formally written as 



= n{rt), the 



= e 



'E^^{rt,z) (FP cavity), 

(15) 

using the approximation Q. The numerical evaluation of ( 15 1 
is the basis what is known as the beam propagation method 
(BPM) ||24)-||26). It has been used for the simulation of both 
passive f27] and active Q, 1 28 1, p9| high-power laser cavities. 

For the case of a spatially and temporal constant index, 
n = const., ( 15 I can be evaluated exactly to yield 



< J G^{r't^rt,z' - z)E^{rt,z)dxdy (FP cavity) (16) 



with 



G±(r;-n,z'-z) 

e(±(z'-z)) 



2tt{z' - z) 



(17) 



where Q denotes the Heaviside step function. The integral 



equation (16 1 together with the propagator (17 1 is known as 
Huygen's integral in the Fresnel approximation. It has been 
successfully applied for the simulation of laser cavities as 
outlined in pO| , | |3T| , including those of semiconductor lasers 

lug-Jig. 

Based on (15 1 and the boundary conditions (|6]l it is possible 
to construct round trip operators M*. One starts at some 
position z = zq within the cavity and performs a full roundtrip. 
Depending whether we start into forward (+) or backward (— ) 
directions, the eigenvalue problems 



M±i;± (rt, zo) = l.mE^{rt, zq) 



(18) 



are obtained with 

Im = e' (19) 

being the eigenvalues. The eigenfunctions of are the mode 
distributions E^{rt,ZQ) at the position zq. 

A very popular method for solving ( fTS] ! is based on the 
Fox-Li approach, cf. |[30| and the references therein. The 
idea is to choose a normalized, more or less arbitrary start 
distribution E^^{rt,zo) and to apply the round trip operator 
recurrently until one arrives (hopefully) at a steady 
state, i.e. a distribution that does not change its pattern in a 
round-trip, except for a reduction in amplitude and a phase 
shift which yields 7,„. Thereby typically the phase factor 
exp( - 2iRe(A^)L - 2i/3oi)) is omitted. 

The Fox-Li approach delivers the mode with the largest 
1 7m I if converged which is, however, not guaranteed. From a 
mathematical point of view, it is related to the power (also 
called vector or von-Mises) iteration to determine the largest 
eigenvalue of a matrix and the corresponding eigenvector 
It is known, that the generated sequences converge only if 
the eigenvalue is simple and well separated from the other 
eigenvalues. Hence, the algorithm fails if there are cavity 
modes having identical or nearly identical dampings lm{ilm). 



(18 1 for passive laser cavities based on Amoldi or Lanczos 
iteration. Recently, an alternative approach based on a finite- 
element method for solving (|9]l directly, without recourse to 
( 18 1, has been proposed |35|. 

Simulations of active laser cavities based on BPMs still 
utilize the Fox-Li approach Q, |36|, which work quite well 
for single-transverse mode lasers. The approach fails for multi- 
transverse mode lasers and at high output power if nonlinear- 
ities such as self-focusing and filamentation start to dominate, 
cf. 1 37 1 for a recent discussion of the stability of the Fox-Li 
approach. The association of the nonconvergence of the Fox- 
Li iteration with an dynamically unstable laser behavior (as in 
p8| ) should be done with care. Instead, for multi-mode high- 
power lasers, a time-dependent approach based on ([T]) should 
be preferred, although numerically more challenging. 



D. Beyond paraxial approximation 

With improving computer capabilities, direct numerical so- 
lutions of Maxwell or Helmholtz equations without recourse 
to the paraxial approximation become more and more fea- 
sible. For example, in |39 | a scalar Helmholtz equation has 
been solved in the {x, z) (lateral, longitudinal) plane for a 
high-power laser, restricted, however, to short cavities. An 
numerical model based on the Helmholtz operator and a first- 
order time derivative obtained from Maxwell equations by 
separating the rapidly varying term exp(iLLii) has been used in 
pOl for the simulation of a monolithically integrated master- 
oscillator power-amplfier. Another possibility is the expansion 
in terms of waveguide modes |41) , |42) considered in the next 
Subsection. 



E. Waveguide modes 

We can expand the left and right traveling waves in terms 
of waveguide modes, 

E^ {rt , z, t) = ^ A± (z, t)xu (r* , z, t) , (20) 



which are solution of 



(21) 



subject to the boundary conditions (|7]l. The waveguide 
modes are exact solutions of Maxwell equations of the form 
E^{rt)exp{iujt — i/3^z) if the dielectric function is translation- 
ally invariant along z ||43). 

Eq. (21 1 is in general a non-hermitian eigenvalue problem 
for the complex propagation factor /3^. An overview of several 
methods to solve (21 1 can be found in p5j. Expansion into 



waveguide modes allows an exact treatment of waveguide 
discontinuities such as Bragg gratings |44| or laser facets |45| 
in order to calculate their reflectivities. In what follows we 



introduce an approximate solution of (21 
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Fig. 1 . Profiles of real part of refractive index(blue solid, left axis), intensity 
(red solid, right axis) and constant phase (red dotted, right axis) of the leaky 
waveguide under study. The real part of the effective index is also shown 
(blue dotted, left axis.) 
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Fig. 2. Modal losses due to radiation into an infinitely thick substrate versus 
the thickness of the cladding layer. The vacuum wavelength is Aq = 1 ^tm. 



F. Effective index method 

For a nearly planar waveguide, i.e. if the mode profile 



time |47|-|49|, it has been recently re-discovered for lasers 



x{x, y) does not vary strongly along the x direction, (21 
be solved perturbationally |46]. We choose 

X{x,y) = ^{x,y)(l){x) 



(22) 



where the dependence of $(a;, y) on x is only parametrically. 
If <i>(a;, y) is the solution of the eigenvalue problem 



kl 



'eff. 



] $ = 



with 



then 



with 



(x) obeys the eigenvalue problem 

dx"^ 



+ [kl {nl, + ^nl,)~P^]4>^Q 



(23) 



(24) 



(25) 



(26) 



In (23 I, n{x,y) is a typically real-valued chosen refractive- 



index distribution of a reference waveguide and neff(a;) is the 
so-called effective index at position x. The correction Arieff 



given by (26 1 could include e.g. modifications of the index due 



to carrier density and temperature effects as well as the imag- 
inary part. The approach sketched above is called "effective- 
index method" in the semiconductor laser community. Eq. (23 1 
is widely used in the design of high-power lasers, because it 
allows the optimization of the layer structure with respect to 
optical confinement, far-field divergence and optical losses. In 
what follows we will consider an example. 

G. Worked example 

Lasers grown on GaAs substrates, suffer from a leakage of 
the lasing mode into the substrate, or a coupling of the lasing 
mode with substrate modes because the refractive index of 
the substrate is typically larger than the effective index of the 
lasing mode. Although this effect is well-known for a long 



grown on GaN substrates |50|. A notably exception are lasers 



with a broad GaAs waveguide core and a multi-quantum well 
active region emitting around 1100 nm |5l]|. 

Assuming an infinite thick substrate, the associated losses 



can be determined by solving ( 23 1 subject to the condition 



d^ 

dy 



'eff 



(27) 



at the boundary between cladding and substrate. Eq. ( [27| ) is 
an example of a boundary condition of type (|8]l. Note that the 
correct branch of the complex square root function has to be 
chosen in order to ensure an outgoing wave. As a consequence 



of (27 1, the effective index rieif becomes complex valued, even 
though there is no absorption in the waveguide. The imaginary 
part of the effective index describes the loss of the modes 
through radiation into the substrate. If there is no absorption in 
the substrate or gain in an active layer, the field diverges with 
increasing distance from the cladding layer If the absorption 
or gain is sufficiently large, this so-called leaky mode becomes 
a proper, albeit complex, guided mode. 

Simulation tools using a PML boundary condition often 
fail to calculate the radiation losses of leaky waveguides 
correctly fST]. Using a transfer-matrix based method, the 
radiation loss can be exactly calculated by solving a complex- 
valued transcendental equation |53|, p4|. We will give here 



the results of a worked example which can be used as a 
benchmark. Fig. [T| shows the profile of the index and the profile 
of the intensity of the mode with the highest effective index 
for a vacuum wavelength of Aq = 1 /im. 

Fig. |2] shows the radiation losses in dependence of the 
thickness of the cladding layer between waveguide core and 
substrate. The dependence is exponential as to be expected, 
except for small values of da, where the mode becomes 
strongly leaky. Due to the fact that the radiation loss increases 
with the effective index, a careful adjustment of the thickness 
of the cladding layer between substrate and waveguide core 
can be used to discriminate higher order vertical modes pT) . 

A second consequence of the optical leakage is the appear- 
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Fig. 3. Vertical profiles of tlie far field intensities for a substrate infinitely 
thick (top) and with a finite thickness (bottom). 



ance of additional peak(s) in the profile of the far field intensity 

2 



Pff{0) oc cos2(e) / ^{yy''°''"'^'^-^^ydy 



(28) 



in dependence on the vertical divergence angle Q. If the mode 
profile in the substrate 



$(2/) oc e-^'=''\/"L-»crf!/ 



is inserted into ( |28| l, there is a resonance if 

sin'(e,) ^n^^b-^eff 



(29) 



(30) 



holds. Depending whether a substrate with an infinite (Fig. 
|3] top) or finite (Fig. |3] bottom) thickness is assumed, one 
resonance or two resonances, respectively, will be observed in 
the far field profile. The magnitude and width of the peaks 
depend on the absorption in the substrate and the bottom 
surface roughness and metalization. The additional peak(s) in 



the farfield can be observed if the condition < n^i^^b — n 



eff < 1 



holds. If rijjjb — n^jf > 1, the internal propagation angle 
Qi = acos(neff/?^sub), which is the angle of the normal of 
the line of constant phase shown in Fig. [T] with respect to the 
z axis, is larger than the critical angle of total reflection and 
no peak appears. 

A third consequence of the coupling of the lasing mode 
with substrate modes is a modulation of the intensity in the 
center of the waveguide core (or the optical confinement 
factor) in dependence on the wavelength as shown in Fig. 
|4] The modulation period depends strongly on the width of 
the waveguide core because the smaller the core the larger the 
wavelength dependence of the index of the mode confined to 
the core. This effect results in a modulation of the modal gain 
and hence in a corresponding modulation of the spectrum of 



the amplified spontaneous emission p9| , |55| and the optical 
spectrum above threshold |56|. 

It should be noted that the GaAs p-contact layer causes 
mode coupling phenomena, too, if the thickness is not properly 
chosen |48|. 




1.00 

wavelength X [nm] 



1.05 



Fig. 4. Relative intensity in the waveguide core versus vacuum wavelength 
for different width of the waveguide core as indicated. 



H. Longitudinal modes 

If we insert (j20]l into (j9]) and neglect the coupling of 
different transverse modes due to a spatially varying group 
index, temporally varying index or spatially varying index 
perturbations, the mode amplitudes obey 



dz 



with the abbreviation 

A/3, = 



[Pi 



Pi] 



For a Fabry-Perot laser, where the coupling coefficients 

^^ klJ^^Eldxdy 
'^'^ 2^0 J E^dxdy 



(31) 



(32) 



(33) 



vanish, Eq. (31 1 can be analytically solved subject to the 
boundary conditions (|6]). From the complex mode frequencies 



and ( 10 1 the wavelengths of the modes relative to the reference 
wavelength can be determined to 



X2 



2LTTn„ 



nk - Ll3o - iRe(A/3,) 

(34) 

where k is the longitudinal mode index and (po and ipL 
are the phases of the reflectivities. The spacing between the 
wavelengths of different transverse modes belonging to the 
same longitudinal mode is approximately given by 



A. 



A. 



Ao 



(35) 



with nmod,i/ = Re(/3,,)/fco. 

Several cases can be distinguished depending on the differ- 
ence of the modal indices governing the wavelength spacing 
of the transversal modes. For example, in a ridge waveguide 
laser with a narrow ridge, the spacing of the wavelengths of 
the lateral modes is in the order of 1 nm which is much 
larger than that of the longitudinal modes. Hence, every lateral 
mode is associated with a comb of longitudinal modes which 
can overlap. In a broad-area laser with a wide ridge or gain 
region, on the other hand, the spacing of the wavelengths of 
the lateral modes is typically much smaller than that of the 
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longitudinal modes, so that in a measured optical spectrum 
each longitudinal mode splits into several peaks belonging to 
the different lateral modes. 

Assuming a ID waveguide with width W and constant real 
refractive index surrounded by perfect electric walls, the 
modal indices are 



with J/ = 1, 2, 3, ■ 



(36) 



Thus the wavelength spacing with respect to the fundamental 
mode I' = 1 is given approximately 



(1 



(37) 



which has been used by several authors fST], f58l for studies 
of broad-area lasers. It should be noted that ( (37] i is an 
approximation, because it neglects the true profile of the index 
and the penetration of the fields into the exterior region. 

III. Time-dependent active cavity simulation 
A. Model equations 

In order to assess the multi-mode behavior, modal instabil- 
ities and filamentation effects of wide-aperture semiconductor 
lasers a time-dependent approach is the most appropriate one 
as outlined in Section II. Due to computer limitations, until 
now only the vertical-projected equations |59|, ||60| 



^'g.eff 



dz 



1 d'^E^ 



dt 



+ 



fcg[r 



'eff " 



'effj 



2/3o 

PI 



2/3o 



E^ + 



2/3o 



-E^ 



spont 



with 



^0 

2/3o 



(38) 



(39) 



have been dealt with. They are obtained by introducing an 
Ansatz like (22 1 into ([T|i and adding a Langevin source i^spom 
on the rhs describing spontaneous emission. Dispersion of the 
imaginary part of the index (optical gain) must be additionally 
included to eliminate the high-fc instability | |6T| , which can 
be done on a microscopic level [62 J , | [63) , as an effective 
polarization p9) modeling a Lorentzian shape of the gain, 
with higher order time derivatives f64'|, by means of a digital 
filter [65 1, or by a convolution integral [66|. 



Eq. ( 38 1 has to be supplemented by equations governing the 
carrier dynamics. In all models so far published this is done on 
the level of a diffusion equation for the excess carrier density 

N, 



dN 

Ik 



J_ 
ed 



(40) 



assuming charge neutrality in the active region. The thick- 
ness of the active region is denoted by d and the rates 
of spontaneous (including non-radiative and radiative) and 
stimulated recombination by R and i?stim, respectively, and 

V = {d/dx,d/dz). 



60 



50 



40 



30 



S 20 



10 



— I — 
P = 7W 



-| — 1 — 1 — 1 — 1 — I — I— 

— w/o thermal lensing 

— with thermal lensing 




lateral angle [degrees] 

Fig. 5. Time-averaged lateral profiles of the far tield intensity of a gain-guided 
BA DFB laser without and with the thermal lensing effect at an output power 
of P = 7 W. 



The current density is given by 167 



j{x,z) 



j [[/ — U'p{x, z)\/r if a;, e active stripe 

[V^C/pj/fi elsewhere 

(41) 

with U being the bias and r and Vl being the resistivity and 
the sheet resistance, respectively, of the p-doped layers. The 
dependence of the current density within the active stripe on 
the Fermi voltage f/p, i-e. the spacing of the electro-chemical 
potentials of electrons and holes, results in a preferred current 
injection into regions with a low carrier density, thus coun- 
teracting spatial holeburning |68J, [69 J . The reason is that 
f/p decreases and thus the difference U — Up increases with 
decreasing N. 

As shown in | |70l , the diffusion coefficient reads 

,d[/p 



D = fip{po + N) 



dN 



(42) 



with /ip being the mobility of the holes in the active layer 
and pq the equilibrium density of the holes. For Boltzmann 
statistics, it is given hy D = 2ksTiJ,p/e (ks Boltzmann 
constant, T temperature, e elementary charge). 
The rate of stimulated recombination is given by 



i?stim = %effgeff + 



with the effective gain 



(43) 



(44) 



The output power Pq l at z ^ and z = L is then given by 



nujoVg,,fid{l-rlr) / \E^{x,z)\l^nLdx. (45) 



Note, that the dispersion of the gain must be included into 
(|43]), too. 

B. Worked example 

The multi-peaked and not diffraction-limited lateral field 
profile of wide-aperture semiconductor lasers has been a long 
standing problem and has been investigated in the past by 



numerous authors |38|, |7I|-|74|. Although the broadening 
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of the far field of continuous-wave operating lasers with 
increasing power can be at least partially attributed to the 
thermal lensing effect pS) 



| |75| , 1 76 1, a complete picture of 
the origin and mechanism has yet to be revealed. 

One mechanism is believed to be due to carrier induced 
antiguiding, i.e. the reduction of the refractive index with 
increasing carrier density. This leads to a self-focusing mech- 
anism because the index increases in regions of high intensity 
due to a depletion of the injected carrier density and can result 



in the formation of lasing filaments 1 7 1 1 . 



As stated in [72 1, the maximum local increase of the 
real part of the effective index is given by Re((5neff) — 
aH5eff,th/2fco because the carrier density can not be depleted 
below transparency. Hence, filamentation effects in state-of- 
the-art high-power strained quantum-well lasers having a low 
alpha factor an around 2 and a small threshold gain .geff,th 
due low internal and outcoupling losses (large cavity length) 
should be of less importance compared to lasers investigated 
earlier 

Another mechanism that could explain the multi-peaked 
structure and the broadening of the farfield is the simultaneous 
lasing of a large number of waveguide modes, originating 
from a built-in or thermally induced waveguide. Indeed, re- 
cent experiments reveal, that even at currents several times 
above threshold the lateral modes can be clearly identified by 
spectrally -resolved near- and farfield measurements p7] , [58) . 

In what follows we present exemplary results obtained with 
the simulation tool "WIAS-LASER" ||60|, ||77), 1^| for a 
broad-area distributed feedback (DFB) laser [79^1 . The device 
has a cavity length of L = 3 mm, a total width of 500 /im 
with an active stripe width of W = 90 and a coupling 
coefficient of k — 3.2 cm^^. 

The discretization steps are Ax = 0.5 /im and Az — 5 /im 
into the lateral and longitudinal directions, respectively. The 
latter corresponds to a time step of At = 31.7 fs. The 
simulation is performed over a total time window of 20 ns, but 
temporal averaging and Fourier transformation is done during 
the last 12 ns thus excluding the turn-on behavior. 



The dependence of the effective gain and index change on 
the carrier density is modeled as ges{N) — g'ln (N/Na-) and 
Aneff(iV) = — V n'N, respectively, with g' — 33 cm^^, A'n- = 
1.4 • 10^^ cm"''', n' = 2.5 • lO"^'^ cm^. The square root like 
dependence of the index on the carrier density is motivated 
by the result of microscopic calculations |80| and was firstly 
used in (28). 

We will consider a device without a built-in index step 
emitting a temporal averaged output power of about P ~ 7 W 
at a current of / = 8 A. Fig. [5] shows the temporal averaged 
profiles of the intensity of the far field 

for two cases, without and with a temperature induced index 
profile due to self-heating. The stationary temperature profile 



has been calculated in advance by the tool JCMsuite 1 76 1, | [81) . 
If the thermal lens is not taken into account, the far field profile 
is asymmetric. This could be caused by a symmetry-breaking 
mechanism inherent to purely gain-guided devices |74|. With 
the thermal-lensing effect the far field has a symmetric shape, 
but is more divergent. 

A deeper insight can be gained from the spectrally resolved 
far fields as shown in Figs. [6] and [7] In Ref. |82| the parabolic 
shaped dispersion curve associated with a longitudinal mode, 
visible in Fig. |6] was misinterpreted as a sign of lateral 
filamentation. We think that the opposite is true: The parabola 
indicates the surviving lateral mode structure, and its blurred 
appearance reflects the filamentation. Due to the lack of any 
carrier-density independent index profile, the modes react very 
sensitively to any changes of the carrier density, caused by 
local intensity fluctuations. 

The index profile resulting from the self-heating stabilizes 
the lateral modes, so that they can be clearly identified in 
the spectrally resolved far field shown in Fig. |7] according to 
their number of nodes and antinodes. The parabolic shaped 



dispersion curve is caused by the fact, that according to ( 35 1 
the lasing wavelengths of the modes decrease and that the far 
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Fig. 8. Time-resolved near field profiles of the gain-guided BA DFB laser 
without the thermal tensing effect at an output power of P = 7 W. 
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Fig. 9. Time-resolved near field profiles of the gain-guided BA DFB laser 
with the thermal tensing effect at an output power of P = 7 W. 



field divergence increases with rising mode index. A similar 
interpretation was also given in |83J . We should note, that 
there is only one longitudinal mode lasing due to the DFB 
operation. 

The time-resolved near fields shown in Figs. [8]and|9]reveal 
the non-stationary behavior, which was already found earlier 
l|8T|. It can be interpreted to be caused by mode beating in 
time and space. Obviously, in the purely gain-guided case 
(Fig. [8]l less modes are involved compared to the case with 
a superimposed index profile (Fig. [9]l. The appearance of the 
bright spots can be thought to be the result of a constructive 
interference of the lateral modes. Although the results are 
encouraging, further analytical and numerical investigations 
are needed in order to reveal the relative contributions of 
lateral mode and filamentary structures to the optical field of 
wide-aperture lasers and to obtain quantitative agreement with 
experimental results. 

IV. Stationary drift-diffusion based simulation 

A. Basic equations 

The standard model |f85l-f87l for the numerical evaluation 
of the distributions of the electron and hole densities N and 
P, temperature T and electron and hole densities j„ and j.^ 
in semiconductor lasers is the energy transport model, which 
consist of the Poisson equation for the electro-static potential 

- V [ssVifi] ^C-P + N (47) 

with es being the static dielectric constant and C the charged 
impurity density, the continuity equations 



Vj„ — R + i?stim 
-Vjp =R + i?stim 

and the heat flow equation 

- V[klVT] = H 



(48) 
(49) 

(50) 



with kl being the thermal conductivity. The various implemen- 
tations of the model differ in the assumptions for the current 



densities the j„, jp and the heat source H. Expressions in 
agreement with the thermodynamic principles can be found in 



Ref. |88|. The rate of stimulated recombination is given here 
by 



(51) 



with Xv obtained by solving ( [2T| ). 

The equations are supplemented by proper boundary condi- 
tions |88|. For example, the boundary condition for the heat 
flow equation reads 



vniyr ^h[T,-T] 



(52) 



where u is the normal unit vector and is the heat sink 
temperature. The heat transfer coefficient h models the heat 
flow to the exterior region. 

The optical output power for a given applied bias U is 
obtained by a solution of (31 1. However, in Fabry-Perot 



lasers it is sufficient to determine the forward and backward 
propagating optical power P+ = I^JP — \A~\'^, 

which are solutions of the equations 



± 



dP^ 



dz 



2Im(/3, 



ao 



(53) 



with Py — P+ + P^ . The loss coefficient ao includes all 
loss mechanisms which are not already accounted for in the 
imaginary part of j3, such as additional internal or scattering 
losses. The equations ( [53] l have to be solved subject to the 
boundary conditions 

P+{{)) = R^p-{0) and P-{L)=RlP+{L). (54) 

The lasing wavelength can be approximately determined 
by determining the maximum of the integral Im(/3y + 
Vty/v^^y)dz. 

Eqs. (31 1 or (53 1 can be conveniently solved by the "Treat 
Power as a Parameter" (TPP) method as introduced in | |68| 
and used in f89l. It is based on the observation, that the 
relative propagation factor in (31 1 or imaginary part 

Im(/3i,) in (53 1 can be considered to be a function of the 
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Fig. 10. Measured (solid red) and simulated optical output power (left axis) 
and conversion efficiency (right axis) versus injection current. The simulations 
were done without (dashed blue) and with (solid green) longitudinal spatial 
holeburning. 



bias, the local power and the unknown wavelengths of the 
lasing and nonlasing modes. Thus, in a first step one calculates 
A/3i, as a function of U, the power of the lasing waveguide 
mode Pjy^ = P+ + P~, and the wavelengths by solving the 
drift-diffusion and waveguide equations in the transverse cross 
section and stores the results in a in a look-up table. In a 



second step, ( 53 i is solved by interpolating in the look-up 



table. For the lasing waveguide mode, lm{n^^ ) has to vanish 
and one has to determine the output power and the wavelength 
for given U. For the non-lasing waveguide modes one has to 
determine il^ for the given U and the power distribution of 
the lasing mode. 

Longitudinal spatial hole burning (LSH) is included auto- 
matically via the power dependence of Im(/3^) in Eq. ([53] l. If 
Im is evaluated at the average power P,y, = / [P+ + 
P^]o?z/L in the cavity, the usual model neglecting LSH 
is recovered, because (53 i is linear and can be analytically 
solved. 

B. Worked example 

In what follows we present results obtained with the sim- 
ulation tool "WIAS-TeSCA" ||90) for a broad-area laser. The 
device with a cavity length of L = 6 mm and an active stripe 
width of T4^ = 95 /xm has a double quantum well (DQW) 
InGaAs/GaAsP active region embedded into 2.2 /xm wide 
Alo.3oGao.7oAs confinement and Alo.85Gao.15As cladding lay- 
ers. The optical confinement factor of the DQW totals 1.8%. 
We performed a one-dimensional simulation in the transverse 



cross section neglecting any lateral effects as outlined in (91 1 



and ||92J. The pre-factor for the gain and the Shockley -Read- 
Hall recombination life times ( t„ = Tp = 1.3 ns identical 
for all layers) have been fitted on the results of length- 
dependent pulsed measurements of threshold current and slope 
efficiency of uncoated devices. The cross-sections of free- 
carrier absorption have been chosen to /cn = 3 x 10^^^ cm^ 




1 2 3 4 5 
longitudinal position z [mm] 



Fig. 1 1 . Longitudinal profiles of optical power (top) and modal gain (bottom) 
at an injection current of / = 25 A without (dashed blue) and with (solid 
green) longitudinal spatial holeburning. 



and fcp — 7 X 10 cm^ for electrons and holes, respectively. 



Fig. 10 shows the measured and simulated power-current 



characteristics of the laser. Experimentally, a maximum output 
power of P = 20 W at a current of / = 25 A was 
achieved limited by thermal rollover. The experimental data 
were reproduced theoretically in a two-stage process with and 
without LSH. The reason for the rollover of the characteristics 
occurring already without LSH has been explained in detail 
in (92]. With increasing current and power, the reduction of 
the gain caused by the temperature rise (AT = 42 K at 
/ = 25 A) must be compensated for by a corresponding 
increase in the carrier densities in the active region, which 
leads as a consequence also to an increase in the carrier 
densities in the confinement layers. Another effect is the 
bending of the quasi-Fermi energy of the holes with increasing 
applied bias and, as a consequence, the bending of the band- 
edge energies due to the voltage drop in the p-doped layers. 
This leads to a corresponding linear increase in the electron 
density with increasing distance from the active region up to 
n w 3 X 10^^ cm^'^. The increased carrier densities in the 
active region and in the bulk layers give rise to enhanced 
non-stimulated recombination and free-carrier absorption 1 ,93) , 
1941. 

However neglecting LSH, not only the roll-over power, but 
also the conversion efficiency is exaggerated which is due 
to the overestimation of the slope-efficiency. Good agreement 
between measurement and simulation is only achieved if LSH 



is included, as Fig. 10 reveals. 

To illustrate the mechanism of LSH, Fig. 11 shows the 



longitudinal profiles of the optical power, the modal gain and 
the injected current density for the two models at a current of 
/ = 25 A. The rear facet with a reflectivity of Rl = 0.98 is 
located at z— 6 mm and the front facet with a reflectivity of 
Po = 0.005 at z = 0. LSH leads to a strong depletion of the 
gain in the vicinity of the front facet, which is compensated by 
a corresponding rise of the gain at the rear facet. This results 
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in a weaker increase of the power towards the front facet and 
hence a lower power at z = 0, i.e. a lower output power, and 
is in contrast to Ref | ,95J where based on an analytical model 
it was found that the values of the power at the facets don't 
dependent on the fact, whether LSH is considered or not. 

It should be noted, that even with LSH there is still a small 
deviation in roll-over power and conversion efficiency between 
measurement and simulation. One reason is the neglect of the 
series resistances of the p-contact and the substrate in the sim- 
ulation, which lead to an additional voltage drop reducing the 
conversion efficiency but also to additional heating reducing 
the roll-over power. 

C. Heterojunctions 

Heterojunctions characterized by discontinuities in the 
edges of the conduction and valance bands as well as abrupt 
changes of effective masses and mobilities require special 
attention. If the validity of the thermodynamic based drift- 
diffusion approach is assumed, the electro-chemical potentials 
should be continuous through the heterojunctions. The rea- 
son is, that the electro-chemical potentials and the inverse 
temperature can be interpreted as Lagrange multipliers in the 
functional for the maximization of the entropy subject to the 
constraints of charge and energy conservation |96| . 

The widely used thermionic emission theory which es- 
tablishes additional conditions for the current densities jOT) , 
| [98| is restricted to situations where the current essentially 
flows perpendicular to the junction and cannot be applied 
to intersections of three or more materials | |99| . In high- 
power lasers there is no need to apply the thermionic emission 
theory with its difficulties because typically all heterojunctions 
(except at the QWs) are graded in order to avoid any drops of 
the electro-static potential. 

D. Quantum effects 

In the simulations presented in subsection IV-B the carriers 
in all layers were treated within the drift-diffasion approach, 
neglecting quantum effects. For more accurate simulations 
electrons and holes confined in active QWs must be treated 
in a special manner As outlined in detail in [87) , |99|- p02) , 
the transport of the confined carriers in the in-plane directions 
{x, z) can be described by classical drift-diffusion, but in the 
perpendicular direction y the carriers have to be described by 
their quantum-mechanical wave functions. Therefore, the car- 
riers population partitions into carriers which have sufficient 
energy to be be considered as unconfined and those which are 
confined. The scattering between both populations is described 
by a capture rate which has to be included into the continuity 
equations ( [48] l as a recombination rate for the unconfined and 
a generation rate for the confined carriers. 

It should be noted that the statistics of both carrier pop- 
ulations is governed by different electro-chemical potentials. 
Whereas the quasi-chemical potentials of the unconfined car- 
riers depend directly on the bias applied to the contacts, 
the quasi-chemical potentials of the unconfined carriers are 
determined by the capture rates. Coupling of the populations 



Although the sketched picture where carriers are described 
quantum-mechanically in one direction and classically in the 
others has been included in several simulation tools, there 
are still open questions which have to be solved. The energy 
which partitions the carrier populations can not be always 
clearly defined, for example in the case of a QW located in 
an unbiased pn-junction. Another issue is the correct treatment 
of the transport in multi-quantum well structures between the 
QWs through the barriers. 

The impact of the non-equilibrium between confined and 
unconfined carriers in high-power lasers on the internal effi- 
ciency, e.g., has been investigated until now using only rate- 
equation approaches | |103| , not the drift-diffusion model. 

V. Summary and Outlook 

We presented models for the calculation of passive cavity 
modes, for the investigation of the spatiotemporal behavior of 
the optical field and of the stationary simulation of the light- 
current characteristics of high-power semiconductor lasers. 
Future work should be directed at full three-dimensional 
calculation of the optical field and an improvement of the 
physical models underlying the time-dependent active cavity 
simulations, including a better description of the carrier and 
heat transport by the energy transport model instead of the 
diffusion equation and stationary heat conduction equation. 
Finally, we would like to mention that the dependence of 
many material parameters and functions such as mobilities 
and absorption coefficients on composition, temperature and 
doping is not well known and needs to be improved. 

Appendix 
Proof of mode orthogonality 

Let us write down ([TJ for 2 modes with indices m and m': 
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(55) 
(56) 

(57) 

■Et, = 



(58) 

Next (ISSll is multiplied with E^^^,, with Et\ with 



£,„, and (58 1 with ~E^. Then all equations are integrated 



over the cavity and added to obtain 

fco iy^m ^m') 



c/3o 



™fig [E+E^, + E^E+,] dxdydz 



takes also place via the Poisson equation (47i 



— [e+e;;^, - e;^e+] dxdydz 

- Vt£;+, - E+^tE;^,] dxdydz = (59) 
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where the derivatives have been factored out. The second and 
third terms vanishes due to the boundary conditions (j6]l-(j7|. 
Hence, the first term has to vanish, too, with leads to pi] ). 
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